pure subroutine pardtc(tx,nx,ty,ny,c,kx,ky,nux,nuy,newc,ier)
! subroutine pardtc takes the knots and coefficients of a bivariate spline, and returns the
! coefficients for a new bivariate spline that evaluates the partial derivative (order nux, nuy) of
! the original spline.
!
! calling sequence:
! call pardtc(tx,nx,ty,ny,c,kx,ky,nux,nuy,newc,ier)
!
! input parameters:
! tx : real array, length nx, which contains the position of the knots in the x-direction.
! nx : integer, giving the total number of knots in the x-direction (hidden)
! ty : real array, length ny, which contains the position of the knots in the y-direction.
! ny : integer, giving the total number of knots in the y-direction (hidden)
! c : real array, length (nx-kx-1)*(ny-ky-1), which contains the b-spline coefficients.
! kx,ky : integer values, giving the degrees of the spline.
! nux : integer values, specifying the order of the partial
! nuy derivative. 0<=nux<kx, 0<=nuy<ky.
!
! output parameters:
! newc : real array containing the coefficients of the derivative.
! the dimension is (nx-nux-kx-1)*(ny-nuy-ky-1).
! ier : integer error flag
!
! restrictions:
! 0 <= nux < kx, 0 <= nuy < kyc
!
! other subroutines required:
! none
!
! references :
! de boor c : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
! dierckx p. : curve and surface fitting with splines, oxford university press, 1993.
!
! based on the subroutine "parder" by Paul Dierckx.
!
! author :
! Cong Ma
! Department of Mathematics and Applied Mathematics, U. of Cape Town
! Cross Campus Road, Rondebosch 7700, Cape Town, South Africa.
! e-mail : cong.ma@uct.ac.za
!
! ..scalar arguments..
integer, intent(in) :: nx,ny,kx,ky,nux,nuy
integer, intent(out) :: ier
! ..array arguments..
real(RKIND), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1))
real(RKIND), intent(out) :: newc((nx-kx-1)*(ny-ky-1))
! ..local scalars..
integer :: i,j,kx1,ky1,lx,ly,l1,l2,m,m0,m1,nkx1,nky1,nxx,nyy,newkx,newky,nc
real(RKIND) ak,fac
! ..
! before starting computations a data check is made. if the input data
! are invalid control is immediately repassed to the calling program.
ier = FITPACK_INPUT_ERROR
if (nux<0 .or. nux>=kx) return
if (nuy<0 .or. nuy>=ky) return
kx1 = kx+1
ky1 = ky+1
nkx1 = nx-kx1
nky1 = ny-ky1
nc = nkx1*nky1
ier = FITPACK_OK
nxx = nkx1
nyy = nky1
newkx = kx
newky = ky
! the partial derivative of order (nux,nuy) of a bivariate spline of degrees kx,ky is a bivariate
! spline of degrees kx-nux,ky-nuy. we calculate the b-spline coefficients of this spline
! that is to say newkx = kx - nux, newky = ky - nuy
newc(:nc) = c(:nc)
if (nux>0) then
lx = 1
x_deriv_order: do j=1,nux
ak = newkx
nxx = nxx-1
l1 = lx
m0 = 1
do i=1,nxx
l1 = l1+1
l2 = l1+newkx
fac = tx(l2)-tx(l1)
if (fac>zero) then
do m=1,nyy
m1 = m0+nyy
newc(m0) = (newc(m1)-newc(m0))*ak/fac
m0 = m0+1
end do
endif
end do
lx = lx+1
newkx = newkx-1
end do x_deriv_order
endif
if (nuy>0) then
ly = 1
y_deriv_order: do j=1,nuy
ak = newky
nyy = nyy-1
l1 = ly
do i=1,nyy
l1 = l1+1
l2 = l1+newky
fac = ty(l2)-ty(l1)
if (fac>zero) then
m0 = i
do m=1,nxx
m1 = m0+1
newc(m0) = (newc(m1)-newc(m0))*ak/fac
m0 = m0+nky1
end do
endif
end do
ly = ly+1
newky = newky-1
end do y_deriv_order
m0 = nyy
m1 = nky1
do m=2,nxx
do i=1,nyy
m0 = m0+1
m1 = m1+1
newc(m0) = newc(m1)
end do
m1 = m1+nuy
end do
endif
return
end subroutine pardtc